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In a previous communication (W.J.T. Bos and J. -P. Bertoglio 2006, Phys. Fluids , 18, 031706), 
a self-consistent Markovian triadic closure was presented. The detailed derivation of this clo- 
sure is given here, relating it to the Direct Interaction Approximation and Quasi-Normal 
types of closure. The time-scale needed to obtain a self-consistent closure for both the en- 
ergy spectrum and the scalar variance spectrum is determined by evaluating the correlation 
between the velocity and an advected displacement vector-field. The relation between this 
latter correlation and the velocity-scalar correlation is stressed, suggesting a simplified model 
of the latter. The resulting closed equations are numerically integrated and results for the 
energy spectrum, scalar fluctuation spectrum and velocity-displacement correlation spectrum 
are presented for low, unity and high values of the Schmidt number. 
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1. Introduction 

The stochastic nature of turbulence introduces the necessity of a statistical descrip- 
tion. Typical statistical descriptors are the mean velocity and the moments of the 
turbulent fluctuations. In isotropic turbulence the mean velocity is zero and the ve- 
locity fluctuations reduce to their geometrically simplest form. Still the interaction 
of a large number of lengthscales renders the description difficult. A phenomeno- 
logical description of isotropic turbulence was proposed by Kolmogorov [Is], leading 
to a prediction of the distribution of energy over different lengthscales in a tur- 
bulent flow assuming that high Reynolds number turbulence is locally isotropic 
and that mode interactions are local in scale-space. Experiments [3], presented at 
the famous Marseille Conference in 1961, confirmed the predicted energy distri- 
bution, proportional to k~^^^, with k the wavenumber. Deriving this description 
from first principles, i.e. starting from the Navier-Stokes equations, is a formidable 
task. Indeed phenomenological models {e.g. references [1,13]) reproduced the cor- 
rect wavenumber dependence of the energy distribution, but the models did not 
have any direct relation to the Navier-Stokes equations. 

A first attempt to build a model for the multi-scale dynamics of a turbulent flow 
could start from the observation that the single-point velocity statistics seem close 
to Gaussian. If is assumed that all velocity moments behave in a Gaussian manner, 
one is directly led to the erroneous conclusion that all different lengthscales are in- 
dependently decaying under the influence of viscosity since the interaction between 
modes is governed by triple correlations which are zero in a Gaussian field. Indeed, 
in a Gaussian field modes are statistically independent and this precludes energy 
transfer among scales. It is thus clear that at the level of third order moments. 
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cumulants, defined as the differences of quantities with respect to their Gaussian 
values, can not be neglected. 

A logical attempt to obtain a correct description of turbulence would then con- 
sist in closing the hierarchy of moments by an assumption of Gaussianity at the 
first level which does allow the triple correlations to be non-zero, in other words, 
at the level of quadruple correlations. This quasi- normal assumption 0-0] is well 
documented in the textbook by Monin and Yaglom @]. Unfortunately, this ap- 
proach leads to non-physical results, since the energy-spectrum predicted by the 
quasi- normal model becomes negative at long-times [9|, as was carefully suggested 
previously by Kraichnan (lol |. 

A great breakthrough in the theory of turbulence was the introduction of the Di- 
rect Interaction Approximation (DIA) pjj]. The great difference with the foregoing 
attempts is that it is a two-time theory, introducing a new quantity, the response 
function, corresponding to the response of a turbulent flow on an infinitesimal 
perturbation. The big difference with the quasi-normal approximation is that cu- 
mulants of all orders are allowed to be non-zero and quadruple moments are not 
assumed to be Gaussian. This issue was at the heart of the discussion between Ian 
Proudman and Robert Kraichnan at the conference in Marseille [i3, [13] (see also 
[3]). Later works showed how these higher order cumulants can be determined 



within the DIA framework [15|, ll6(]. The DIA theory suffered however from both a 



theoretical and a practical weakness. 

The theoretical weakness consists in the fact that DIA does not predict the 
correct inertial range behavior (i.e. in agreement with Kolmogorov's prediction). 
This was traced back to the fact that its dynamical behavior is not invariant to 
Galilean transformations 17|] and the defect was corrected for by reformulating the 
theory in a Lagrangian coordinate system [l^. This yielded a description, called 
the Lagrangian History DIA which gives correct asymptotic behavior for the energy 
spectrum. However, the formulation becomes extremely complicated in this new set 
of coordinates. A simplified, or abridged, version was proposed in the same work 
and an even simpler descri ptio n of the same type is the Lagrangian Renormalized 
Approximation (LRA) [li,]!^]. 

The practical weakness of these theories, both DIA and the Lagrangian versions, 
is the two-time character which precludes a computation of long-time statistics, 
unless the time-history is truncated at some time difference. To obtain a practi- 
cal statistical model for the dynamics of turbulence, one needs a procedure called 
Markovianization. This procedure is based on the assumption that the statistical 
two-time behavior of the turbulent quantities is known, so that the theories can 
be simplified to a dynamical description which only depends on the current time. 
Kraichnan proposed such a procedure and applied it to DIA, yielding an approx- 
imation known as the Test Field Model ^]. The two-time quantities are in this 
approach assumed to be exponentially decaying with a typical time-scale which is 
determined by comparing the decorrelation of a fully compressible field, the test 
field, with the decorrelation of the advecting velocity field. A simplified version of 
the Test Field Model is the Eddy- Damped Quasi-Normal Markovian model [13], 
in which this time-scale is heuristically modeled by a phenomenological straining 
time-scale of the Heisenberg type [3]. 

In a previous work fi^] we proposed a method to determine this time-scale self- 
consistently in the framework of the EDQNM theory, using the observation that 
this time-scale in Lagrangian History DIA is the correlation-time of a fiuid particle 
along a trajectory. In the present article we will show how our theory is related 
to the DIA. This will in particular show why a certain number of terms which 
appeared in the original formulation have to be taken equal to zero for consistency 
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with DIA. We will also apply the approach to the problem of the diffusion of a 
passive scalar. In the framework of the EDQNM theory, we will propose a model 
for the scalar flux spectrum, which significantly simplifies the level of complexity 
with respect to the models proposed in previous works lii-i^]. 



2. Lagrangian dynamics, Markovianization and determination of the 
correlation time-scale 

Let us introduce the spectral tensor, which is the key quantity in the statistical 
description of a turbulent flow. 

$ij(fc,t, S) = T(r) ^Ui{x,t)uj{x + r, s)| (1) 
= Ui{k,t)u*{k,s) (2) 

with a Fourier transform with respect to r , t) the Fourier transform of the 
velocity and the overline denoting an ensemble average. The conjugate u*-{k^s) = 
Uj{—k, s) since the velocity is a real quantity in physical space. In incompressible, 
mirror-symmetric, isotropic turbulence this quantity is entirely determined by a 
scalar function, E{k,t,s), the energy spectrum, by the relation 

^,^ik,t,s) = ^^^^P.j{k), (3) 

with 

P^,{k)=6.,-^. (4) 

The goal of the present work is to obtain a statistical model for the dynamics of 
E{k, t, t) = E{k, t). The original DIA theory was based on the above quantities, or 
representatives, yielding a k~^l'^ inertial range scaling. This behavior was related to 
the fact that an Eulerian description does not distinguish between the decorrelation 
of the velocity modes by the effects of pressure and viscosity on one hand and 
the sweeping of the small scales by the large scales on the other hand. Sweeping 
does not contribute to the energy transfer among scales which is the physical 
mechanism behind the k~^/^ prediction. Indeed an Eulerian observer would have 
no way to see if the decorrelation is due to sweeping or not. To correctly capture the 
dynamics one needs a Lagrangian description. In other words, one needs to change 
the representative from ([T|) to a Lagrangian formulation. Kraichnan reformulated 
his theory in Lagrangian coordinates [IS]. Hereto he introduced the generalized 
velocity, Ui{x,t\s), defined as the velocity measured at time s on the trajectory 
which passes at t through x. The first time argument is called the labeling time, 
and the second, the measuring time. This velocity is illustrated in Figure [H We can 
associate the Lagrangian position function to this velocity, Xi{x^ t|s), defined as the 
position of the fluid particle at time s on the trajectory which passes through x at 
time t. These quantities are related to the Eulerian velocity and coordinate system 
by, 

u{x, t\t) = u{x, t), (5) 
x{x, t\t) = X. (6) 
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The key quantity that we wih use in the following is the Eulerian displacement 
field X related to the generalized velocity field through the relation 

Xi{x,t) = / Ui{x,t\s)ds = Xi{x,t\t) — Xi{x,t\0). (7) 

^0 

The displacement field is an Eulerian field, which measures at each point x in a 
fixed, Eulerian, coordinate system the displacement with respect to its original 
position (at t = 0) of the particle which is at point x at time t (see for example the 
discussion in reference [23]). The displacement vector is also illustrated in Figure 
[H Its evolution is given by 

'-^^ + Uj{x,t)—Xi{x,t)=Ui{x,t). (8) 

Note that since the displacement field is a quantity measured in an Eulerian coor- 
dinate system, the evolution equation contains a convective derivative. 

Using the generalized velocity, the spectral velocity-tensor now contains four 
time- arguments, 



^ij{k,t\s;t'\s') = Ui{k,t\s)u*{k,t'\s'). (9) 

However, in the Lagrangian Renormalized Approximation [19] or in the Abridged 
Lagrangian History DIA 18|] , the resulting model equations only contain quantities 
of the form 



<^ij{k,t\t;t\s) = Ui{k,t\t)u*{k,t\s), (10) 

i.e., one of the velocities is evaluated at a fixed point, since the labeling and mea- 
suring time coincide, and the other is moving along its Lagrangian trajectory. The 
resulting model thus only contains two-time quantities, like the Eulerian DIA, but 
the time-difference is evaluated along a particle trajectory. The difference between 
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LRA and the Abridged Lagrangian History DIA is that in the former, time deriva- 
tives are taken with respect to the measuring time, whereas in the the latter these 
derivatives are taken with respect to the labehng time. 

We note that the Lagrangian velocity is not solenoidal. The Eulerian velocity 
is incompressible so that ^ij{k,t\t\t\s) is solenoidal in its first index i. Isotropy 
and mirror-symmetry imply then that the tensor ^ij{k,t\t]t\s) is proportional to 
Pij{k). We can therefore define the spectrum, 

t X E(k, t\t;t\s) „ , , , 

^i.j{k,t\t;t\s) = ^^J^; ' ^ Pi,{k), (fl) 

but this relation is not valid for ^ij{k,t\s;t'\s'). Since only two-time quantities are 
used in the following, the spectral tensors and spectra will not be denoted by the 
four time-arguments, but will be denoted by their short-hand notation, 

^ij{k,t\s) = <^>ijik,t\t;t\s) (12) 
E{k,t\s) = E{k,t\t;t\s). (13) 



This notation should not introduce confusion, because we will only consider spectra 
in which three of the time- arguments are evaluated at t. Note that the Lagrangian 
and Eulerian representative coincide when s = t, 



E{k,t\t) = E{k,t,t) = E{k,t). 



(14) 



The Abridged Lagrangian History DIA constitutes a closed set of equations, yield- 
ing a inertial range behavior. The equation governing the evolution of E{k, t) 
is 



d ,9 
ot 



E{k,t)= ff [ [xy + z^]pE{q,t\s)[k'^G{k,t\s)E{p,t\s) 

J J A Jo 

dpdq 



-p'^G{p,t\s)E{k,t\s)\ ds- 



pq 



(15) 



in which G{k, t\s) is the response function measured along a trajectory. The symbol 
A indicates the domain in the pq-plane in which k,p,q can form a triangle (in other 
words \p — q\ < k < \p + q\) and x, y, z are given by 

X = -Piqi/{pq) 
y = kiqi/{kq) 

z = kiPi/{kp). (16) 



The LRA equations differ from (|15p in that the RHS contains quantities of the form 
G{k,s\t) and E{k^s\t) instead of quantities of the form G{k,t\s) and E{k,t\s). In 
other words, the labeling and measuring times are inverted in the two-time quan- 
tities. We will not consider the LRA in the following. In the Abridged Lagrangian 
History DIA, equation ([15]) is part of a system of three equations, the other two giv- 
ing the evolution for E[k,t\s) and G{k,t\s)^ respectively. The equations, although 
simpler than the original Lagrangian History DIA equations, remain rather com- 
plicated and depend on the entire time-history of the flow. Our goal is to simplify 
the description to obtain a model that does not depend on the entire time-history 
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of the flow, while retaining as much as possible from the original description. We 
will now explain how we obtain such a model. 
We introduce the correlation function R{k,t\s) as, 



R{k,t\s) 



E{k,t\s) 
E{k,t,t)' 



which allows to rewrite (jlSh as 
d 



dt 



E{k,t) 



[xy + z^]pE{q,t) @{kpq)k'^ E{p,t) 



with 



Q{kpq) 



-e{jpkq)p^G{p,t)E{k,t) 



G{k,t\s)R{p,t\s)R{q,t\s)ds 



dpdq 
pq 



(17) 



(18) 



(19) 



At this point we have not introduced any assumptions with respect to the Abridged 
Lagrangian DIA. A first assumption that we make is 



R{k,t\s) = G{k,t\s) for t > s. 



(20) 



which is a fluctuation-dissipation hypothesis, which leads to 



E{k)= jj Q{kpq)[xy + z^]pE{q) [eE{p) - p'^E{k)] (21) 
with 

@{kpq)= [ G{k,t\s)G{p,t\s)G{q,t\s)ds. (22) 

The fluctuation-dissipation theorem can be proven to hold rigorously in equilib- 
rium thermodynamics, but in turbulence, corresponding to an out-of-equilibrium 
dissipative system, this is not the case [1^. It is observed in [sO] that in the 
Abridged Lagrangian History DIA, the fluctuation-dissipation relation is violated. 
In the Lagrangian Renormalized Approximation, however, which is different from 
Lagrangian History DIA with respect to the choice of the representative, the 
fluctuation-dissipation relation holds surprisingly [l^, [13] . The precise form of the 
relation between G and R is therefore not uniquely determined and can have im- 
portant consequences for the Markovian models which result from the assumption, 
in particular when wave- like phenomena are considered [sH . 

In the above description, only the triad time-scale Q{kpq) depends on the time- 
history of the flow, through the temporal dependence of the response function 
G{k,t\s). Therefore, by modeling the time-dependence of G{k,t\s), or equivalently 
R{k,t\s), a single-time description can be obtained. It is at this point that we will 
proceed in simplifying the description, by introducing specific assumptions on the 
time-dependence of R{k,t\s). A plausible assumption is that two-time correlations 
decay exponentially with the time-difference. This cannot be correct for very small 
times since the time-dependence should be differentiable at zero time difference and 
such a discontinuity can influence the results in the case of compressible turbulence 
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[33] • A Gaussian time-dependence could be assumed to get correct behavior around 
zero. It was however shown that in the incompressible case the results are relatively 
insensitive to the exact choice (33| , since only the integrals over the exponential 
are required and the exact form of the time-correlation does not appear in the 
final expressions. For convenience we will assume therefore an exponential time- 
dependence here: 



R{k, t\s) = exp 



\t - s\ 



[ r{k,t)\ 



(23) 



Note that also the choice that r is here a function of t (and not of s) is a nontrivial 
assumption. We can integrate both sides of this expression with respect to s to 
obtain an expression for T{k,t), 



T{k,t) 



R{k,t\s)ds I 1 — exp 



|t-to| 

T{k,t) 



(24) 



For the sake of simplicity, we will consider here the case in which the turbulence 
has been created a time long before we evaluate it, by taking the limit of to tending 
towards —00. This gives the expression 



r(/c,t) 



R{k, t\s)ds 



E{k,t\s] 
-00 E{k,t\t) 



ds. 



(25) 



The essential step to get a closed model expression for T{k, t), which was proposed 
in [2^, was the explicit integration of E{k,t\s) with respect to s. Using relation 
d?]), we can evaluate the integral (1251) . 

/E{k,t\s)ds = 27rk'^J'(^r) \ / Ui{x,t\t)ui{x + r,t\s)ds 
-00 w— 00 

= 2TTk'^T(^r) ^Ui{x, t\t)Xi{x + r, f)| 



Fik,t,t), 



(26) 



in which we introduced a new quantity, the velocity-displacement cross-correlation 
spectrum F{k, t, s), 



F{k,t,s) 



and we will write 



Ui{k,t)Xj{-k,s) = ^^^2 -^^i(^) 



F{k,t,t) = F{k,t). 



(27) 



(28) 



At this point we have obtained an Eulerian single-time expression for the La- 
grangian time-scale. 



rik,t) 



F{k,t) 
E{k,t)' 



(29) 



If the equation for E{k, t) is completed by an equation for F{k, t), we obtain a closed 
set of equations, if the equation for F{k, t) does not introduce other unknown 
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quantities. In |23|] this was accomplished within the framework of the EDQNM 
closure. In the next section we will obtain similar expressions, using the DIA for- 
malism which will clarify uncertainties in the original derivation. It was observed 
by Kraichnan f[§\ that the evolution ol E(k,t\t) in the Abridged Lagrangian His- 
tory DIA is given by the same equation as E{k, t) in the Eulerian DIA. The only 
difference is that the two-time quantities are now evaluated along particle trajecto- 
ries. We will use this as an ansatz in the following and we will derive the Eulerian 
DIA expression for F(k,t), in which we will evaluate the two-time correlations 
using Kraichnan's generalized Lagrangian description. In section [6] the resulting 
equations will be numerically integrated and results will be presented for both 
the dynamics of the kinetic energy spectrum and the dynamics of a passive scalar 
advected by turbulent flow. 



3. Derivation of the evolution equation for the velocity-displacement 
cross-correlation 

We will focus in this section on the derivation of the evolution equation for the 
wavenumber spectrum of the correlation between the velocity m and the displace- 
ment vector Xi of a fluid particle. The DIA formalism will be used to evaluate the 
triple correlations which appear in this equation. A Markovianization procedure 
is applied to the resulting expressions to obtain a single-time expression. We will 
try to give a self-contained description of all the operations used in this derivation. 
However, for a justification of the assumptions in the DIA we refer to the original 



work 111] or to later works that discuss the justification of the Lagrangian vari- 
ants of the DIA The statistically isotropic and mirror-symmetric case 
will be considered but in principal this derivation could also be performed for the 
anisotropic case. 
The starting point are the equations for Ui and Xi, 

^l^i^ + u-(x t)—u (x t)--^^^ + u—u {x t) 

dJU^ ^ dX^ ^ ^^.(^,i)^X,(a:,t) = u,{x,t), (30) 

which in Fourier space are written as 

+ uk'^]ui{k,t) = -^Pijm{k) J j 6k-p-qUj{p,t)um{q,t)dpdq 
dXi{-k,t) 



dt 



with 



kj I i 5k~p-qUj{-p,t)Xi{-q,t)dpdq + Ui{-k,t), (31) 



Pijm{k) — kjPim{k) -\- k„iPij[k). (32) 



Multiplying the first equation by Xi(—k,t), the second by Ui{k,t), summing and 
ensemble averaging gives 

+ uk^]F{k, t, t) = T'^ik, t) + T^{k, t) + E{k, t, t) (33) 
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with 

T'^{k,t) = -mk^Pijmik) 11 5k^p-qUj{p,t)um{q,t)Xi{-k,t)dpdq 



T^{k,t) = 2mk^kj j J 5k-p-qUi{k,t)uj{-p,t)Xi{-q,t)dpdq. (34) 

The DIA is a perturbation method in which the small parameter is the influence of a 
single triad interaction on the velocity field, compared to the velocity field obtained 
by the interaction of all other triad interactions. In practice, DIA expressions are 
obtained by replacing, in the correlations that one is interested in, the velocity field 

by 



Ui ^ ^ + ) (35) 



and using a similar expression for the displacement vector Xi. The part uf'^ is the 
velocity field that would exist if all triad interactions are retained, except for one 
that we will call {k,p, q). Since in a homogeneous fiow an infinite number of triads 

interact, uf^ is supposed to be close to the actual velocity field. The contribution 

u"P corresponds to the velocity field induced by the interaction {k,p, q) only. The 
justification of this approach can be interpreted in terms of the coupling strength 
between the difi'erent modes that should be sufficiently small. If this coupling is 
weak enough and the number of modes tends to infinity, the DIA procedure should 
give good results. This was tested in reference [33] on a toy-model in which both 
the coupling strength and the number of modes could be varied. 
After replacing all u's and X's in by ()35p and developing, we find that 



uf\p)u?J{q)xl'\-k) + 



zeroth order 



u^^\p)u^^\q)X^''\-k) + uf\p)ug\q)X^^\-k) + uf\p)u^^\q)X^'\-k) 



first order 

+ second and higher order] dpdq, (36) 

in which time-arguments were omitted to shorten the expression. A similar expres- 
sion is obtained for . In the DIA only the first non- vanishing order is retained 
in addition to the zeroth order contribution. Higher order contributions should 
contribute insignificantly if the contribution of a single triad is indeed infinitesi- 
mally small. For the triple correlations considered here, the zeroth order contri- 
bution (first term on the RHS) in the above expression is equal to zero. The first 
nonvanishing contribution is the first order contribution. The Direct Interaction 
contributions can be written for, for example Uj{p,t) and Xi(k,t), as 



r-t 

p,t) = -iPjabip) G{p,t,s) u'^^\k,s)u^^\-q,s 

Jo L 

) i-k, t) = zka f (fc, t, s) [4°) (-P, s)xf^ i-q, s) + (-p, s)] ds, 

Jo L J 

(37) 
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in which G{k,t,s) and {k,t, s) are the response function of the velocity and 
displacement field, respectively. These are unknown quantities at this point. 

These direct interaction contributions are substituted into expression (j36p . To 
illustrate the derivation of the final expression we take the first term on the second 
line of expression (j36p . which we will denote Ti(k), and we will work out the final 
expression for this term. Upon substituting ()37p in this term we obtain, 



(fc, i-q, s)u'^^ {q, t)xf^ {-k, t)dpdqds. (38) 

This expression is exactly what would be obtained using the quasi-normal ap- 
proach, for a particular choice of G, the viscous Green-function. However, DIA does 
not only take into account viscous decorrelation, but also non-linear decorrelation 
and the DIA expression for G will therefore differ from the viscous Green-function. 
In the limit, which we will consider, in which the flow is constituted of a very 
large number of modes, we can write the correlation in ([38|) as a product of double 
correlations (note that this is no assumption of Gaussianity of the velocity modes, 
as for example discussed in [37|), 



(39) 

in which we used the isotropic expressions for the velocity and velocity- 
displacement tensors. Introducing this in the above expression we find 

Ti(A:) = -^k^ mi 6k-p-qG{p, t, s)Pi,m{k)Pjab{p)Pbm{q)Pai{k) ^^/^y ^ ^^^^'j^2 dpdqds, 

(40) 

which can be written as 

Ti{k) = -^k^ JJ ^k-p-qj ds G{p,t,s)kp fi{k,p,q)E{q,t,s)F{k,s,t){16Tr^k\'^)-^dpdq 

(41) 

with fi{kpq) = {kpy^ Pajm{k)Pjab{p)Pbm{q.) in which we used that 
Pijm{k)Pai{k) = Pajm{k). The Convolution integral can be rewritten as a func- 
tion of the wavenumber only. The details of this simplification are given in the 



appendix of the book by Leslie [38(]. The resulting expression is 

j j 6k-p-qdpdq —7- J J 2'Kpqk~^ dpdq (42) 

in which A indicates the domain in the pg-plane in which /c,p, q can form a triangle 
(in other words |p — (?| < k < \p + q\)- This yields 

Tiik) = -\l f dsG{p,t,s)h{k,p,q)p''E{q,t,s)F{k,s,t)^^. (43) 
8 7a Jo P Q 
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We see that the unknown quantities in this description are the 2-time quantities 
E{k,t, s), F{k, s,t),G{k,t, s). It is here that we replace these quantities by their 
analogues defined on Lagrangian trajectories, 

E{k,t,s) E{k,t\s) 
G{k,t,s) G{k,t\s) 
F{k,t,s) F'{k,t\s) 
F(k,s,t) ^ F"(k,t\s), (44) 

with 

F'{k,t\s) = 27^ k'^Ui{k,t)Xi{-k,t\s) 

F"{k,t\s) = 2-Kk^Ui{k,t\s)Xi{-k,t) (45) 

in which the difference between F' and F" arises due to the asymmetric character 
of the velocity-scalar displacement correlation. As in section [21 expression (f23l) . we 
will use the assumption that the two-time quantities decay exponentially in time 
by posing 



G{k, t\s) = H{t - s) exp[-(t - s)7]{k)] (46) 

G^{k,t\s) = H{t - s) exp[-{t - s)r]^ (k)] (47) 

E{k, t\s) = E{k, t)[G{k, t\s) + G{k, s\t)] (48) 

F'{k, t\s) = F{k, t)[G{k, s\t) + G{k, t\s)] (49) 

F"{k,t\s) = F{k,t)[G^{k,s\t) + G^{k,t\s)]. (50) 



The Heaviside functions H(t) appear since the response functions correspond to 
the response at time s to a perturbation applied at time t. There is obviously 
no response before the perturbation is applied. Note further that (l48]) . (l49]) and 
(150p . which link the response functions to the velocity time-correlations, are, like 
expression ()20p . a specific form of the fluctuation-dissipation theorem. Using the 
above expressions, the time-integral will, in the long-time limit, yield 



t 

Ft 



ds G{p,t\s)E{q,t\s)F"{k,t\s) = {kpq)E{q,t)F{k,t) 

(51) 

with 



&^{kpq) = (52) 

V^ik) +r]{p) +r]{q) 



We repeat the foregoing analysis for the remaining two terms in and the three 
terms of . The Direct Interaction contribution stemming from the Xi in each 
term yields two contributions, whereas every Ui yields only one contribution (c./. 
the symmetry used in expression (00])). The total number of terms is then eight. 
In the following all quantities are evaluated at time t so that we will omit the time 
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arguments. The resulting equation is 

[^^+uk']F{k) = Y,T,ik) + Eik) (53) 

i=l 

with 



Ti (k) = 


1 

~8 


/ Q^(kva)fT(k V a)v^E(a)F(k^ 
J A 


dp dq 
p q 






1 

~8 


/ Q^(kr)a)f-,(k v a)a^ E(v)F(k'' 
J A 


dp dq 
p q 




T-^(k) = 


1 

^8 


/ Q^(kna)U(k v a)k^E(n)F(a^ 
J A 


dp dq 
p q 




Ti{k) = 


1 

+i 


r 

/ @^{kpq)fi{k,p,q)k^E{q)F{p] 
J A 


dp dq 
p q 


(57) 


nik) = 


~4 


[ e^{qpk)f5{k,p,q)p^E{k)F{q] 
J A 


dp dq 
p q 


(58) 


nik) = 


+i 


[ e^{qpk)Mk,p,q)k^E{p)F{q] 
J A 


dp dq 
p q 


(59) 


Tj{k) = 


~4 


[ Q''iqpk)Mk,p,q)q^E{k)F{p] 
J A 


dp dq 
p q 


(60) 


Tsik) = 


~4 


[ e^{qpk)fs{k,p,q)q^E{p)F{k] 
J A 


dp dq 
p q ■ 


(61) 



Symmetries allow to show that Ti{k) = ?2(A;) and T^i^k) = T/^i^k), which reduces 
the number of terms to six. The geometrical functions /i, .., /s can be expressed as 
a function of the cosines x,y,z, defined in expressions (fT6]l . This yields 



/l 


= {kp)-'Pi23{k)P2iA{p)P3A{q) = 


2{xy + z"^] 




(62) 


/2 


= {kq)-^Pi23{k)P24{p)P3iAiq) = 


2{xz + y'^) 




(63) 


/3 


= {k)-^hP2u{k)Pn{p)P2A{q) = 


1 — xyz — 




(64) 


U 


= {k)-^kiP234{k)P23{p)Pu{q) = 


1 — xyz — 


2y^z^ -y^ + z^ 


(65) 


h 


= {kpy^kiP23{k)PMp)PM{q) -- 


= xy + z^ - 


- z + zx^ 


(66) 


k 


= {k)-^kiP2M{k)Pi3{p)P2A{q) = 


h 




(67) 




= {kq)-^kiP2z{k)Pn{p)q2 = z{x 


+ yz) 




(68) 


fs 


= 2{kqy^kiPi2{p)q2 = 2{y + xz] 






(69) 



Note that fi and /2 are not needed when symmetry is used to show that Ti{k) = 
T2{k) and T3{k) = Ti{k). Since /g = /s, we eventually only need five different 
factors. 

The last step in obtaining closure is the determination of ri{k) and 'q^{k). Com- 
parison of equation ()48p with ([23]) shows that within the framework of our assump- 
tions, 

ri{k) = T{k)-^ = E{k,t)/F{k,t). (70) 
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Note that the current model, as it was presented in (23| defined r/(/c) as 

ri{k) = E{k,t)/F{k,t) +vk^. (71) 

Indeed it is now clear that within the framework of the assumptions introduced in 
the derivation of the curent model, the viscous frequency should not be added. It 
is implicitly contained in the expression for F{k). 

We have obtained a closed expression, if the unknown quantity rj^{k) is specified. 
We note here that the displacement is a quantity which will not diffuse along a 
trajectory, neither by pressure, nor by molecular effects. This suggests the choice 

r?^(A;) = 0, (72) 



in analogy with Lagrangian History DIA arguments for the passive scalar. There 
is however a difference between the equation of the displacement vector and the 
passive scalar. Indeed, the equation of the displacement vector ([30|) contains on the 
RHS the velocity. This forcing is correlated with the displacement. Some analogy 
exists between this case and the linearly forced Burgers equation. In the latter 
case the triad relaxation time is directly determined by the correlation-time of the 
velocity, as was shown by applying the LHDIA technique to Burgers' equation [3^ . 
In order to test the influence of this time-scale on the results, we will also perform 
simulations in which 

i^^{k) = T-\ (73) 

in which T is the integral timescale of the velocity field. It will be shown that this 
choice does not significantly change the results. 



4. Equations for the evolution of the energy spectrum and the scalar 
variance spectrum 

Using the above method, equations can be derived for the evolutions of the energy 
spectrum and of the scalar variance spectrum. The latter is defined by 



Ee{k, t, s) = 2Trk^9{k, t)e{-k, s). (74) 

The resulting equations are identical to the ones given by the EDQNM approach 
or the test-field model. The only difference is the specification of the time-scale. 
The equations are 



at 
at 



E{k) = jj e{kpq)[xy + z^]pE{q) [^Eip) - p'^E(k)] (75) 
Eg{k)= [[ Q\kpq)[l - z']pEg{q) [k^E{p) - p'E{k)] (76) 

JJa pq 



The time-scales which appear in these expressions are 



Using the arguments from the last section, the timescale ri{k) ^ is given by expres- 
sion (|29p. The time-scale associated to the decorrelation of a passive scalar on a 
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Lagrangian trajectory can be inferred from the equation for the scalar, 

(if) 

^ = aA«. (78) 

which suggests that the scalar decorrelation is given by 

r]\k) = ak'^, (79) 

since on a Lagrangian trajectory the scalar only decorrelates under the action of 
diffusion. We note that this choice for the scalar is inherent to the choice of the 
representative ^-^]. For example, the decorrelation time-scale for the velocity is 
modified if the DIA is applied to the strain-rate instead of to the velocity. However, 
we will use the above choice, expression ([79l) . 

For comparison, it is perhaps useful here to recall the expressions for the time- 
scales which are used in the EDQNM model ji^]. These are 



ij{k) = xJj^ s^E{s)ds + uk'^, r]'{k) = X'Jj^ s'^E{s)ds + uk^ 



r]^{k) = Xe^ s^E{s)ds + ak^. (80) 

The value of A was determined and given the value 0.36 by comparison with the 
test-field model, the latter being calibrated by comparison with the DIA in a case 
in which DIA was assumed to give correct results. The value for Ag was given 
the value in (43| . using similar arguments as above regarding the Lagrangian 
decorrelation of a scalar fluctuation. In the same work A' was given a value different 
from A. The parameter A' appearing in Q^{kpq) was determined by comparison 
with experimental results for the turbulent Prandtl number and was given the 
value unity. The different values for A and A' in the two different time-scales are 
in contradiction with the outcome of DIA. However, a certain flexibility is allowed 
if one drops the fluctuation-dissipation theorem. Even more freedom is obtained if 
one changes the representative used in the analysis. From a practical point of view 
one might argue that the appearance of adjustable constants in a model allows the 
user to choose them in comparison with observation rather than by using physical 
arguments and assumptions. 

These expressions for the time-scales used in the EDQNM closure were only 
recalled here for comparison. They will not be used in the remainder of this work. 



5. EDQNM expression for the spectrum of the velocity-scalar 
cross-correlation 

It is at this point that we want to point out a possible simplification of the EDQNM 



model as proposed and derived in references [244271]. In these works the EDQNM 
expression was derived for the spectrum of the velocity-scalar correlation spectrum. 
The resulting expressions contained a time-scale of the form (j80p . However in the 
light of the present derivation it seems that the level of sophistication of that model 
was higher than the standard EDQNM model. Indeed, the considerations in the 
present work suggests a much simpler expression to model the scalar flux spectrum. 
Using the analogy between a passive scalar generated by a mean-scalar gradient 
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(a) 



10" 10' 10"'^ 10"' lO'* 10'^ 10"^ 10"' 10° 10' 



(b) 



Figure 2. Energy spectrum (a) and in compensated form (b) for decaying isotropic turbulence, modeled 
by the present closure approach. 





(a) 



(b) 



Figure 3. Velocity-displacement cross-correlation spectrum (a) and in compensated form (b) for decaying 
isotropic turbulence, modeled by the present closure approach. 



and the displacement- vector, expression (|29]) suggests the fohowing model for the 
scalar flux spectrum Fuoik), 



Fue{k) = TE{k)hF{k) 



TE{k) 



^fJ J s'^ E{s)ds + {v + a)k^ 



-1 



(81) 



"with r the strength of the mean scalar gradient and Xp a model constant. This 
modeling "would be of the level of approximation of the EDQNM model and "would 
avoid having to program the lengthy evolution equations "which appear in the pre- 
vious "works. 



6. Numerical results 



In this section "we numerically integrate equations ()53p . ()75p and (j76p . In these 
expressions we "will insert for r]{k) and r]^{k): 



7]{k) = E{k,t)/F{k,t); rj^ {k) = 0; rf{k) = ak 



(82) 
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(a) (b) 

Figure 4. Energy transfer (a) and spectrum of the typical (eddy-damping) time-scale compared to a 
Heisenbcrg type time-scale, usually employed in the EDQNM model (b). 




(a) (b) 

Figure 5. (a) Passive scalar variance spectrum corresponding to a decaying passive scalar at unity Schmidt 
number. The velocity field is characterized in Figure [2] (b) Compensated form. 



We will first consider decaying turbulence, starting from an initial distribution 

E{k) = Sfc^e-^C^A^)' (83) 

with B chosen to normalize the energy to unity, ki = 10 and the first wavenum- 
ber of the domain is A;o = 1. Wavenumber space is discretized using a geometrical 
discretization ki = /cor*"^. This possibility in two-point closures to use a geomet- 
rical discretization allows to attain much higher Reynolds numbers than in direct 
numerical simulations for a much lower cost. The discretization in the present com- 
putations is approximately 15 wavenumbers per decade. The viscosity is equal to 
u = 8 - 10~®. From the initial condition we let the turbulence decay until it reaches 
a self-similar state in which e/e ~ with the kinetic energy and viscous dissipation 
respectively determined by 



e = J E{k)dk, ^ = / k^E{k)dk. (84) 
All results are evaluated in this phase and the Reynolds number, defined as 
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(c) (d) 

Figure 6. Passive scalar variance spectrum corresponding to a stationary passive scalar, injected in the 
large scales at small and high Schmidt numbers. Scalar spectrum at Sc = 0.01 (a) and in compensated 
form (b). Sc = 20 (c) and in compensated form (d). 



is at the time of evaluation 5000. The scalar spectra are initialized to have the same 
shape as the kinetic energy. Mixing at three different Schmidt numbers, Sc = u/a, 
is evaluated, Sc = 0.01,1,20. For the computations at Sc = 0.01 and Sc = 20, 
in order to clearly produce scaling behavior, both the velocity and the scalar field 
are forced at the lowest wavenumbers, k < 2. For the case in which 5c = 20, The 
Reynolds number is taken equal to R\ = 100 and the discretization is refined to 
50 wavenumbers per decade to more acurately capture nonlocal interactions which 
play an essential role in the physics of the Batchelor range. 



6.1. Results for the velocity field for high Reynolds number decaying 
turbulence 

Results for decaying isotropic turbulence at Rx = 5000 are shown in Figure [2j We 
observe that the Kolmogorov constant is of order 1.7. It is observed that in the 
present model the bottleneck effect, that is, the increasing value of the compensated 
plot in the beginning of the dissipation range, is absent. Note that this bottleneck is 
far more pronounced in the EDQNM model [4J]. Indeed, the bottleneck is generated 
by the fact that nonlocal interactions are cut off in the dissipation range, thereby 
decreasing the transfer of energy so that it piles up . This effect is very sensitive 
to the relative value of the viscous time-scale compared to the local (in /c-space) 
turbulent time-scale. For example, the bottleneck can be increased by decreasing 
the value of A in the EDQNM model. 
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The constant in 

F{k) = Cpe^/^k-^/^ (86) 

is order Cp ~ 2.5 as shown in Figure [3j The energy transfer spectrum is shown 
in Figured] (a). In Figured] (b) we see that if we compare the T]{k) given by our 

model to X^J s'^E{s)ds, there is a short range in which the two are approximately 
proportional, with a value of A ~ 0.5 . 



6.2. Results for the mixing of an isotropic passive scalar 

In Figure [5] we show the spectrum of a passive scalar which is freely decaying. The 
initial spectrum is equal to the spectrum of the velocity field and the spectra are 
computed at the same time- instant as in Figure [2] A clear inertial-convective range 
is observed in which the scalar variance distribution is proportional to k~^^^, in 



agreement with classical arguments [45], l46(]. In the compensated representation it 



is observed that the dimensionless constant in 

Ee{k) = Ceeee-i/3A:-5/3 (37) 

is around 0.3 according to our closure. This is about one half of the generally 
observed value in experiments and numerical simulations ^] . This underprediction 
is a well-known deficiency of Lagrangian History DIA and its relatives. It can be 
cured for by changing the representative on which the closure is based US-Si. 



This would constitute an interesting perspective, but is not tried in the present 
work. We note that EDQNM predicts a more realistic value for Cg We stress, 
however, that this value is directly related to the parameter A'. If this value was 
chosen equal to A in the cited reference, then the value for Cg would be of the 
same order as the one predicted by our theory. We further observe that a large 
bottleneck is observed in all these results. 



6.3. Mixing of an isotropic passive scalar at small and large Schmidt 
number 

We also carried out computations for small and large Schmidt number and the 
results are shown in Figure [6] It is observed in Fig. [6] (a) and (b) that the closure 
manages to reproduce the inertial diffusive range [48,] . 

Eg{k) = \c,ege'l^a-^k-'y\ (88) 

In Fig. [6] (c) and (d), we observe that also that the viscous-convective range is 
reproduced [49|, 

Egik) = C'gegu^/h-^/^k-\ (89) 



6.4. Influence of the correlation time-scale of the displacement vector 

To test the influence of the correlation time-scale of the displacement vector on 
the results, we performed a simulation, in which we changed the time-scale of the 
displacement vector from zero (rj^ik) = 0) to the integral time-scale {rj^ (k) = 
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(a) (b) 

Figure 7. Influence of the correlation time-scale of the displacement vector on the energy spectrum. Energy 
spectrum (a). Black corresponds to r)-^ (k) = 0. The red graph corresponds to rj-^ {k) = T~^. Compensated 
spectrum (b). 



T~^), as described in the last paragraph of section [3l The results are shown in 
Figure \7\ Clearly, the influence of this change in correlation time is not dramatic 
in the current application. Only in the large scales of the energy spectrum a small 
influence is observed. In the inertial and dissipation range the change in time-scale 
does not seem to influence the results. The reason for this is that in Q^{kpq), only 
one of the three wave-vectors, k is influenced by this change. The r}{p) and r]{q) 
rapidly dominate over the r]^{k), since r]{k) is an increasing function of k in the 
inertial range. 



7. Concluding remarks 

In the present work we derived a self-consistent Markovian closure starting from the 
DIA formalism and using the fluctuation-dissipation assumption. Closure was ob- 
tained by relating the turbulence time-scale to the spectrum of the cross-correlation 
of fluid-particle-displacement and velocity, F{k,t). The general idea was announced 
in a letter [23:1], but the details of the idea were not given there. In the present 
work we show how the present closure is related to the Abridged La- 
grangian History DIA, and how the introduction of the displacement 
vector allows to measure the Lagrangian decorrelation in an Eulerian 
reference frame. We want to stress that such a single-time closure could 
perhaps also be obtained applying a similar procedure to the two-time 
correlations in the Lagrangian Renormalized Approximation, which was 
not tried here. 

The closure was numerically integrated at large Reynolds number (Rx = 5000). 
Results were presented for the energy spectrum and transfer spectrum. Results for 
the mixing of an isotropic passive scalar were shown for Schmidt numbers 0.01, 1 
and 20 respectively and classical scaling was reproduced. 

In starting from DIA it becomes clear why a certain number of terms should 
be neglected in previous works. These terms were present in work on the velocity- 



scalar cross-correlation spectrum j24l427l| and correspond in these studies to non- 
linear transfer terms which were linearly dependent on the mean-scalar gradient. 
They should be removed if consistency is required with first-order renormalized 
perturbation theories such as DIA. They appear as terms which are second order 
in the formal expansion parameter in the DIA approximation. 



Furthermore, a simplified model for the scalar flux spectrum was proposed ((ST 
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with a level of complexity comparable to the standard EDQNM model for the 
energy spectrmn. 

Let us finish by noting the following: during the Marseille congress 50 years 
ago, Kraichnan presented his DIA [13] to the scientific community. Many things 
have changed since then: computational power has increased enormously, making 
possible the direct numerical simulation of laboratory experiments, advanced ex- 
perimental techniques allow to visualize almost every feature of a turbulent flow. 
However, the theoretical description of turbulence, derived from the Navier-Stokes 
equations, has still not advanced enormously beyond DIA related approaches. The 
physical insights obtained from these triadic closures are still an indispensable 
building block of our description of turbulent flows. 
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